tic
mcvar=zeros(51,1);
for ii=1:10
    randseed=randn(5000,5000);
    samplepaths=oupath(randseed,1,1,1,1);
    samplepaths=samplepaths.^2;
    mcvar=mcvar+sum(midpoint(samplepaths).^2,2);
end
mcvar=mcvar/50000;
toc

% a=1;b=1;m=1;yic=1;
% time=(0:50)*0.5/50;
% exactvar=b^2/2/a*(1-exp(-2*a*time))'+(m+(yic-m)*exp(-a*time)).^2';
plot(time,mcvar)
% hold on
% plot(time,exactvar,'r')